*Replace serif as Helvetica instead of "Times"
graph set eps fontfaceserif "Helvetica"
graph set eps

use "${clean_data}\gfs_cleaned_merged_database.dta", clear
cd "${output}\"

**Figure 1--Simple bivariate tables
drop good_dad good_mom

gen good_dad=4 if father_relatn==1
replace good_dad=3 if father_relatn==2
replace good_dad=2 if father_relatn==3
replace good_dad=1 if father_relatn==4
gen good_mom=4 if mother_relatn==1
replace good_mom=3 if mother_relatn==2
replace good_mom=2 if mother_relatn==3
replace good_mom=1 if mother_relatn==4
gen good_parent=good_dad +good_mom
tab good_parent

sum thriving [aw=weight]

sum parent_died  [aw=weight]
sum parent_died if gdp_pc_max>=25000 [aw=weight]
sum parent_died if gdp_pc_max<25000 [aw=weight]

cor PCRQ REincome_feelings tertiary hhinc if gdp_pc_max>=20000 [aw=weight]
cor PCRQ REincome_feelings tertiary hhinc if gdp_pc_max<20000 [aw=weight]

sum PCRQ if parent_died==1  [aw=weight]
sum PCRQ if parent_died==0  [aw=weight]
sum PCRQ  if gdp_pc_max>=20000 [aw=weight]
sum PCRQ if gdp_pc_max<20000 [aw=weight]

sum wb_fiveyrs wb_today thriving

collapse (semean) SEthriving=thriving (mean) thriving  [aw=weight], by(good_parent)
foreach v in thriving  {
gen UP=`v'+(SE`v'*1.96)
gen LO=`v'-(SE`v'*1.96)
}

label define label_pq 2 "Very bad relationship with both parents" 3 " " 4 "Parent average is somewhat bad" 5 " " 6 "Parent average is somewhat good" 7 " " 8 "Very good relationship with both parents",  add
label values good_parent label_pq
levelsof good_parent, local(pqval)

levelsof good_parent, local(pqval)
twoway (dot thriving good_parent, horizontal mcolor(blue*.7) msymbol(Dh) msize(large) ndots(0) dotextend(no)  ///
ylabel(`pqval', valuelabels angle(0) labsize(small)) ///
lwidth(thick) color(blue*.7) fintensity(0)) ///
|| rcap UP LO  good_parent, lpattern(dash) lwidth(thick) lcolor(red*.7) horizontal  ///
legend(pos(6)) legend(row(1) order(1 "Thriving in wellbeing" 2 "95% CI") ) ///
|| scatter good_parent thriving, msymbol(none) mlabp(10) mlabel(thriving) mlabsize(small) mlabcolor(black) mlabformat(%9.2f) ///
xlabel(, labsize(medium)) ///
ytitle("") xtitle("Percent thriving in wellbeing", size(small)) ////
title("A. Thriving in wellbeing and parent-child relationship qualtiy", span size(medium)) ///
saving(fig1_dot, replace) ///
plotregion(color(white)) graphregion(color(white))  

twoway (dot thriving good_parent, horizontal mcolor(blue*.7) msymbol(Dh) msize(large) ndots(0) dotextend(no)  ///
ylabel(`pqval', valuelabels angle(0) labsize(small)) ///
lwidth(thick) color(blue*.7) fintensity(0)) ///
|| rcap UP LO  good_parent, lpattern(dash) lwidth(thick) lcolor(red*.7) horizontal  ///
legend(pos(6)) legend(row(1) order(1 "Thriving in wellbeing" 2 "95% CI") ) ///
|| scatter good_parent thriving, msymbol(none) mlabp(10) mlabel(thriving) mlabsize(small) mlabcolor(black) mlabformat(%9.2f) ///
xlabel(, labsize(medium)) ///
xlabel(0(.1).5) xmtick(0(.1).5)  ///
ytitle("") xtitle("Percent thriving in wellbeing", size(small)) ////
title("A. Thriving in wellbeing and parent-child relationship qualtiy", span size(medium)) ///
saving(fig1_dot0, replace) ///
plotregion(color(white)) graphregion(color(white))  

*Figure 2
use "${clean_data}\gfs_cleaned_merged_database.dta", clear

foreach x in  svcs_father svcs_mother {
gen REV`x'=4 if `x'==1
replace REV`x'=3 if `x'==2
replace REV`x'=2 if `x'==3
replace REV`x'=1 if `x'==4
}
gen CUMsvcs = REVsvcs_father + REVsvcs_mother
tab CUMsvcs 

gen vgood_parents=1 if father_relatn==1 & mother_relatn==1
replace vgood_parents=0 if father_relatn>1 & mother_relatn>1
replace vgood_parents=. if father_relatn>4 | mother_relatn>4

sum  vgood_parents CUMsvcs [aw=weight]

collapse (semean) SEPCRQ=PCRQ SEvgood_parents=vgood_parents (mean) PCRQ vgood_parents [aw=weight], by(CUMsvcs)
foreach v in PCRQ vgood_parents {
gen UP`v'=`v'+(SE`v'*1.96)
gen LO`v'=`v'-(SE`v'*1.96)
}

label define label_pr 2 "Neither parented attended religious services" 3 " " 4 "Parents average <1/month" 5 " " 6 "Parents average 1-3 times/month" 7 " " 8 "Both parents attended at least once per week",  add
label values CUMsvcs label_pr
levelsof CUMsvcs, local(prval)
drop if CUMsvcs==.


twoway (dot vgood_parents CUMsvcs, horizontal  mcolor(blue*.7) msymbol(Dh) msize(large) ndots(0) dotextend(no)  /// 
ylabel(`prval', valuelabels angle(0) labsize(small)) ///
lwidth(thick) color(blue*.7) fintensity(0)) ///
|| rcap UPvgood_parents LOvgood_parents  CUMsvcs, lpattern(dash) lwidth(thick) horizontal  ///
legend(pos(6)) legend(row(1) order(1 "Very good relationship with parents" 2 "95% CI") ) ///
|| scatter CUMsvcs vgood_parents, msymbol(none) mlabp(1) mlabel(vgood_parents) mlabsize(small) mlabcolor(black) mlabformat(%9.2f) ///
xlabel(, labsize(medium)) ///
xlabel(0(.2)1) xmtick(0(.2)1)  ///
ytitle("") xtitle("Percent with very good relationship with parents", size(small)) ////
title("B. Parent-child relationship qualtiy and parent religious attendance", span size(medium)) ///
saving(fig2_dot, replace) ///
plotregion(color(white)) graphregion(color(white))  
graph export "Fig_Good_Parents_and_Religiosity.png", replace width(900) height(900)

*Dot originate at zero
twoway (dot vgood_parents CUMsvcs, horizontal   mcolor(blue*.7) msymbol(Dh) msize(large) ndots(0) dotextend(no)  /// 
ylabel(`prval', valuelabels angle(0) labsize(small)) ///
lwidth(thick) color(blue*.7) fintensity(0)), yscale(range(0 1)) xscale(range(0 1)) ///
|| rcap UPvgood_parents LOvgood_parents  CUMsvcs, lpattern(dash) lwidth(thick) horizontal  ///
legend(pos(6)) legend(row(1) order(1 "Very good relationship with parents" 2 "95% CI") ) ///
|| scatter CUMsvcs vgood_parents, msymbol(none) mlabp(1) mlabel(vgood_parents) mlabsize(small) mlabcolor(black) mlabformat(%9.2f) ///
xlabel(0(.2)1) xmtick(0(.2)1)  ///
xlabel(, labsize(medium)) ///
ytitle("") xtitle("Percent with very good relationship with parents", size(small)) ////
title("B. Parent-child relationship qualtiy and parent religious attendance", span size(medium)) ///
saving(fig2_dot0, replace) ///
plotregion(color(white)) graphregion(color(white))  

**FIGURE 1
graph combine  fig1_dot.gph fig2_dot.gph , ///
imargin(small) row(2) col(1) iscale(*.75) ///
plotregion(color(white)) graphregion(color(white)) 
graph export "Fig1_2_Parents_Wellbeing_and_Religion.png", replace width(1800) height(1200) 

graph combine  fig1_dot0.gph fig2_dot0.gph , ///
imargin(small) row(2) col(1) iscale(*.75) ///
plotregion(color(white)) graphregion(color(white)) 
graph export "Fig1_2_Parents_Wellbeing_and_Religion_Origin0.png", replace width(1800) height(1200) 

************
**Multilevel Models
*************
use "${clean_data}\gfs_cleaned_merged_database.dta", clear
cd "${output}\"

global country_level_controls "v2x_polyarchy zwvs_traditional zwps_index lngdp"
global country_valid "zmort_rate"
global demo "male other_gender secondary tertiary hhinc  age_group2 age_group3 age_group4 age_group5 age_group6 age_group7 foreign_born  employed never_married married live_alone no_children one_child two_children  rural village suburb"
global parent_demo "parents_married living_structure2 living_structure3 living_structure4  child_poverty was_abused"
global extended_controls "parent_relig_index self_relig_index REincome_feelings"

foreach x in $country_level_controls {
gen X_`x'=`x'*parent_relig_index
gen REL_`x'=`x'*PCRQ
}

global interaction_religion "X_v2x_polyarchy X_zwvs_traditional X_zwps_index X_lngdp"
global interaction_pcrq "REL_v2x_polyarchy REL_zwvs_traditional REL_zwps_index REL_lngdp"
gen Xparent_relig_index_PCRQ=parent_relig_index*PCRQ

label var Xparent_relig_index_PCRQ "PCRQ X parental religiosity index"

label var X_v2x_polyarchy "Interaction-parental religiosity X electoral democracy"
label var  X_zwvs_traditional "Interaction-parental religiosity X secular values"
label var  X_zwps_index  "Interaction-parental religiosity X Women's rights"
label var  X_lngdp "Interaction-parental religiosity X GDP per capita"

label var REL_v2x_polyarchy "PCRQ X electoral democracy"
label var  REL_zwvs_traditional "PCRQ X secular values"
label var  REL_zwps_index  "PCRQ X Women's rights"
label var  REL_lngdp "Interact-PCRQ X GDP per capita"

***********
*Finding 1: Predict flourishing--GDP and tradition only country level controls
**********

eststo clear

foreach Y in HOPE HAPPY  {
eststo: quietly mixed `Y'  PCRQ $demo  [pw=weight] || iso:, mle nolog vce(cluster iso)  
estadd local log_likelihood `e(ll)', replace
eststo: quietly mixed `Y'  PCRQ $demo $parent_demo  [pw=weight] || iso:, mle nolog vce(cluster iso)  
estadd local log_likelihood `e(ll)', replace

eststo: quietly mixed `Y'  PCRQ $demo $parent_demo lngdp $extended_controls [pw=weight] || iso:, mle nolog vce(cluster iso)  
estadd local log_likelihood `e(ll)', replace
eststo: quietly mixed `Y'  PCRQ $demo $parent_demo lngdp $extended_controls REL_lngdp [pw=weight] || iso:, mle nolog vce(cluster iso)  
estadd local log_likelihood `e(ll)', replace

eststo: quietly mixed `Y'  PCRQ $demo $parent_demo zwvs_traditional $extended_controls  [pw=weight] || iso:, mle nolog vce(cluster iso)  
estadd local log_likelihood `e(ll)', replace
eststo: quietly mixed `Y'  PCRQ $demo $parent_demo zwvs_traditional $extended_controls REL_zwvs_traditional [pw=weight] || iso:, mle nolog vce(cluster iso)  
estadd local log_likelihood `e(ll)', replace

eststo: quietly mixed `Y'  PCRQ $demo $parent_demo v2x_polyarchy $extended_controls  [pw=weight] || iso:, mle nolog vce(cluster iso)  
estadd local log_likelihood `e(ll)', replace
eststo: quietly mixed `Y'  PCRQ $demo $parent_demo v2x_polyarchy $extended_controls REL_v2x_polyarchy [pw=weight] || iso:, mle nolog vce(cluster iso)  
estadd local log_likelihood `e(ll)', replace


eststo: quietly mixed `Y'  PCRQ $demo $parent_demo zwps_index $extended_controls  [pw=weight] || iso:, mle nolog vce(cluster iso)  
estadd local log_likelihood `e(ll)', replace
eststo: quietly mixed `Y'  PCRQ $demo $parent_demo zwps_index $extended_controls REL_zwps_index [pw=weight] || iso:, mle nolog vce(cluster iso)  
estadd local log_likelihood `e(ll)', replace

}

esttab using "T1_predict_wellbeing_multi_level_model.csv", ///
stats(log_likelihood N k, label("Log likelihood" "Sample Size" "k number of parameters")) ///
b(%15.3f) se(%15.3f) label  star(* 0.05 ** 0.01 *** 0.001)   nogaps replace

*re-do with complete statistics
esttab using "T1_predict_wellbeing_multi_level_model_t_ci_p.csv", ///
stats(log_likelihood N k, label("Log likelihood" "Sample Size" "k number of parameters")) ///
cells(b(fmt(3)) t(fmt(3) par) ci_l(fmt(3))&ci_u(fmt(3)) p(fmt(4) par))  ///
nogaps   nolines  collabels(none)  star(* 0.05 ** 0.01 *** 0.001) label numbers replace ///
addnotes("rows correspond to beta values; t-statistics; 95% confidence intervals; p-values.")

***********
*Finding 2: Predict PCRQ--GDP and tradition only country level controls
**********

eststo clear

foreach Y in PCRQ  {
eststo: quietly mixed `Y' $demo  [pw=weight] || iso:, mle nolog vce(cluster iso)  
estadd local log_likelihood `e(ll)', replace
eststo: quietly mixed `Y' $demo $parent_demo  [pw=weight] || iso:, mle nolog vce(cluster iso)  
estadd local log_likelihood `e(ll)', replace
eststo: quietly mixed `Y' $demo $parent_demo parent_relig_index [pw=weight] || iso:, mle nolog vce(cluster iso)  
estadd local log_likelihood `e(ll)', replace

eststo: quietly mixed `Y'  $demo $parent_demo lngdp $extended_controls [pw=weight] || iso:, mle nolog vce(cluster iso)  
estadd local log_likelihood `e(ll)', replace
eststo: quietly mixed `Y'  $demo $parent_demo lngdp $extended_controls X_lngdp [pw=weight] || iso:, mle nolog vce(cluster iso)  
estadd local log_likelihood `e(ll)', replace

eststo: quietly mixed `Y'  $demo $parent_demo zwvs_traditional $extended_controls [pw=weight] || iso:, mle nolog vce(cluster iso)  
estadd local log_likelihood `e(ll)', replace
eststo: quietly mixed `Y'  $demo $parent_demo zwvs_traditional $extended_controls X_zwvs_traditional [pw=weight] || iso:, mle nolog vce(cluster iso)  
estadd local log_likelihood `e(ll)', replace


eststo: quietly mixed `Y'  $demo $parent_demo v2x_polyarchy $extended_controls [pw=weight] || iso:, mle nolog vce(cluster iso)  
estadd local log_likelihood `e(ll)', replace
eststo: quietly mixed `Y'  $demo $parent_demo v2x_polyarchy $extended_controls X_v2x_polyarchy [pw=weight] || iso:, mle nolog vce(cluster iso)  
estadd local log_likelihood `e(ll)', replace


eststo: quietly mixed `Y'  $demo $parent_demo zwps_index $extended_controls [pw=weight] || iso:, mle nolog vce(cluster iso)  
estadd local log_likelihood `e(ll)', replace
eststo: quietly mixed `Y'  $demo $parent_demo zwps_index $extended_controls X_zwps_index [pw=weight] || iso:, mle nolog vce(cluster iso)  
estadd local log_likelihood `e(ll)', replace


}

esttab using "T2_predict_relationship_multi_level_model.csv", ///
stats(log_likelihood N k, label("Log likelihood" "Sample Size" "k number of parameters")) ///
b(%15.3f) se(%15.3f) label  star(* 0.05 ** 0.01 *** 0.001)   nogaps replace

*re-do with complete statistics
esttab using "T2_predict_relationship_multi_level_model_t_ci_p.csv", ///
stats(log_likelihood N k, label("Log likelihood" "Sample Size" "k number of parameters")) ///
cells(b(fmt(3)) t(fmt(3) par) ci_l(fmt(3))&ci_u(fmt(3)) p(fmt(4) par))  ///
nogaps   nolines  collabels(none)  star(* 0.05 ** 0.01 *** 0.001) label numbers replace ///
addnotes("rows correspond to beta values; t-statistics; 95% confidence intervals; p-values.")



***********
*Findings 3: OLS--parent effect by country
*******
use "${clean_data}\gfs_cleaned_merged_database.dta", clear
cd "${output}\"

/*
reg HOPE PCRQ $demo living_structure2 living_structure3 living_structure4 [aw=weight] if country_name=="South Africa"

replace selfid1=. if selfid1<0
bysort country_name: egen max_race=max(selfid1)
gen any_race_var=1 if max_race>0 & max_race!=.
replace any_race_var=0 if max_race==0 | max_race==.
sum any_race_var
tab country_num if any_race_var==0
**5,10,17,18
*/

encode country_name, gen(country_num)

gen B1=.
gen B2=.
gen P1=.
gen P2=.
gen SE1=.
gen SE2=.
gen DF1=.
gen DF2=.

forval i=1/22 {
reg HOPE PCRQ $demo living_structure2 living_structure3 living_structure4 [aw=weight] if country_num==`i'
replace B1=_b[PCRQ] if country_num==`i'
replace SE1=_se[PCRQ] if country_num==`i'
replace P1=_r_p[PCRQ] if country_num==`i'
replace DF1=_r_df[PCRQ] if country_num==`i'
}

forval i=1/22 {
reg PCRQ $demo parent_relig_index [aw=weight] if country_num==`i'
replace B2=_b[parent_relig_index] if country_num==`i'
replace SE2=_se[parent_relig_index] if country_num==`i'
replace P2=_r_p[parent_relig_index] if country_num==`i'
replace DF2=_r_df[parent_relig_index] if country_num==`i'
}

format P1 P2 %15.4f
format DF1 DF2 %15.0f

collapse (count) obs=PCRQ (first) country_name B1 B2 SE1 SE2 P1 P2 DF1 DF2 iso country_num v2x_polyarchy zwvs_traditional zwps_index lngdp zmort_rate ///
wps_education wps_employment wps_fin_inclusion wps_cell_phone wps_parliament wps_no_legal_discrimination wps_justice wps_mortality_ratio wps_son_bias wps_partner_violence wps_safety wps_terrorism wps_conflict ///
(mean) PCRQ thriving HOPE HAPPY ACCEPTED_BY_GOD FINANCIAL parent_relig_index self_relig_index parents_married parent_died  flove_no_dad mlove_no_mom child_poverty was_abused [aw=weight], by(country)

gen T1=B1/SE1
gen T2=B2/SE2

sort country_name
save "country_level_effect_estimates.dta", replace

***********
*FIGURES--Rank country-effects with 95% CIs
************
use "country_level_effect_estimates.dta", clear

gen UP1=B1+(1.96*SE1)
gen LO1=B1-(1.96*SE1)
gen UP2=B2+(1.96*SE1)
gen LO2=B2-(1.96*SE1)

drop if B1==.
drop if country==.
gsort -B2
edit
gsort -B1
edit

gen gdp=exp(lngdp)
sort country
order country	obs	PCRQ	HOPE	HAPPY	ACCEPTED_BY_GOD	FINANCIAL	thriving parent_relig_index	self_relig_index  gdp v2x_polyarchy zwvs_traditional zwps_index B1	B2	SE1	SE2	T1	T2	UP1	LO1	UP2	LO2
export delimited country	obs	PCRQ	HOPE	HAPPY	ACCEPTED_BY_GOD	FINANCIAL	thriving parent_relig_index	self_relig_index gdp v2x_polyarchy zwvs_traditional zwps_index B1	B2	SE1	SE2	T1	T2	UP1	LO1	UP2	LO2	///
using "G:\Survey_Methodology\People\Family Study\Global Flourishing Study\generated data\Country_Level_Effects.csv", replace

egen rank1=rank(B1), field
egen rank2=rank(B2), field

* -labmask- is from the Stata Journal 
labmask rank1, values(country_name)
labmask rank2, values(country_name)

/*
tab rank
twoway scatter rank B1, mcolor(blue*.7) msymbol(D) msize(large) legend(off) yla(1/22, ang(h) grid val noticks) ///
|| rcap UP LO rank, mcolor(green*.7)  msymbol(d) msize(small) horizontal ytitle("") ysc(reverse) ///
plotregion(color(white)) graphregion(color(white)) xtitle("Effect of Parent-child Relationship on Flourishing")

twoway scatter rank B1, mcolor(blue*.7) msymbol(D) msize(large) legend(pos(6)) legend(row(1) order(1 "Mean effect" 2 "95% CI") ) ///
xsc(alt) yla(1/22, ang(h) grid val noticks) ///
|| rcap UP LO rank, mcolor(green*.7)  msymbol(d) msize(small) horizontal ytitle("") ysc(reverse) ///
plotregion(color(white)) graphregion(color(white)) xtitle("Effect of Parent-child Relationship on Flourishing")
*/

twoway scatter rank1 B1, mcolor(blue*.7) msymbol(D) msize(large) legend(pos(6)) legend(row(1) order(1 "Mean" 2 "95% CI") ) ///
 yla(1/22, ang(h) grid val noticks) ///
|| rcap UP1 LO1 rank1, mcolor(green*.7)  msymbol(d) msize(small) horizontal ytitle("") ysc(reverse) ///
saving(c_level_effect1, replace) ///
plotregion(color(white)) graphregion(color(white)) title("A. Effect of Parent-child Relationship on Flourishing", span size(med))
graph export "Country_Level_Effect_Size_PCRQ_on_Wellbeing.png"

twoway scatter rank2 B2, mcolor(blue*.7) msymbol(D) msize(large) legend(pos(6)) legend(row(1) order(1 "Mean" 2 "95% CI") ) ///
  yla(1/22, ang(h) grid val noticks) ///
|| rcap UP2 LO2 rank2, mcolor(green*.7)  msymbol(d) msize(small) horizontal ytitle("") ysc(reverse) ///
saving(c_level_effect2, replace) ///
plotregion(color(white)) graphregion(color(white)) title("B. Effect of Parent Religiosity on Parent-child Relationship", span size(med))
graph export "Country_Level_Effect_Size_Religion_on_PCRQ.png"

**FIGURE 2
*grc1leg is an external package; please install
grc1leg c_level_effect1.gph c_level_effect2.gph   ///
, imargin(2 2 2 2) row(3) col(2) legendfrom(c_level_effect1.gph) iscale(*.75) plotregion(color(white)) graphregion(color(white))  
graph export "FIG_Two_Country_Level_Effects_Plot.png", replace width(2000) height(1300) 

*******
*FIGURES--**Correlation plots***
********
use "${clean_data}\gfs_cleaned_merged_database.dta", clear
merge 1:1 country_name using "${output}\religious_diversity_by_country.dta"

cor B1 B2 PCRQ wps_education wps_employment wps_fin_inclusion wps_cell_phone wps_parliament wps_no_legal_discrimination wps_justice wps_mortality_ratio wps_son_bias wps_partner_violence wps_safety wps_terrorism wps_conflict
cor B1 PCRQ  v2x_polyarchy zwvs_traditional zwps_index lngdp zmort_rate thriving HOPE HAPPY ACCEPTED_BY_GOD FINANCIAL parent_relig_index self_relig_index parents_married parent_died  flove_no_dad mlove_no_mom child_poverty was_abused
cor B2 PCRQ  v2x_polyarchy zwvs_traditional zwps_index lngdp zmort_rate thriving HOPE HAPPY ACCEPTED_BY_GOD FINANCIAL parent_relig_index self_relig_index parents_married parent_died  flove_no_dad mlove_no_mom child_poverty was_abused
cor DIV B1 B2 PCRQ self_relig_index zwvs_traditional

*Test if religious diversity matters as much as religiosity
gen inter_div_rel=self_relig_index*DIV
reg B1 zwvs_traditional DIV
reg B2 zwvs_traditional DIV
reg B1 DIV
reg B2 DIV
reg B1 zwvs_traditional 
reg B2 zwvs_traditional 
reg B1 self_relig_index DIV 
reg B2 self_relig_index DIV 

label var parent_relig_index "Parental religiosity index"

corr PCRQ  parent_relig_index
local corr : di %05.2g r(rho) 
twoway lfit ( PCRQ  parent_relig_index) || scatter ( PCRQ  parent_relig_index), mlabel(iso) ///
ytitle("Parent-child relationship quality") plotregion(color(white)) graphregion(color(white))  ///
saving(relate1, replace) ///
note(r =`corr') title("A. Parent-child relationship and parental religiosity") legend(pos(6)order(1 "Linear Fit" 2 "Observed" )  region(style(none))  rows(1) span) 

label var lngdp "GDP per capita (PPP-adjusted, natural log)"
corr PCRQ  lngdp 
local corr : di %05.2g r(rho) 
twoway lfit ( PCRQ  lngdp) || scatter ( PCRQ  lngdp), mlabel(iso) ///
ytitle("Parent-child relationship quality") plotregion(color(white)) graphregion(color(white))  ///
saving(relate2, replace) ///
note(r =`corr') title("B. Parent-child relationship and GDP per capita") legend(pos(6)order(1 "Linear Fit" 2 "Observed" )  region(style(none))  rows(1) span) 

corr B1  parent_relig_index
local corr : di %05.2g r(rho) 
twoway lfit ( B1  parent_relig_index) || scatter (B1  parent_relig_index), mlabel(iso) ///
ytitle("Relationship effect on wellbeing") plotregion(color(white)) graphregion(color(white))  ///
saving(relate3, replace) ///
note(r =`corr') title("C. Relationship effect on wellbeing and parental religiosity") legend(pos(6)order(1 "Linear Fit" 2 "Observed" )  region(style(none))  rows(1) span) 

label var lngdp "GDP per capita (PPP-adjusted, natural log)"
corr B1  lngdp 
local corr : di %05.2g r(rho) 
twoway lfit (B1  lngdp) || scatter (B1  lngdp), mlabel(iso) ///
ytitle("Relationship effect on wellbeing") plotregion(color(white)) graphregion(color(white))  ///
saving(relate4, replace) ///
note(r =`corr') title("D. Relationship effect on wellbeing and GDP per capita") legend(pos(6)order(1 "Linear Fit" 2 "Observed" )  region(style(none))  rows(1) span) 

**predict religiosity effect

label var zwvs_traditional "Inglehart secular values index"
corr B2  zwvs_traditional
local corr : di %05.2g r(rho) 
twoway lfit ( B2  zwvs_traditional) || scatter (B2  zwvs_traditional), mlabel(iso) ///
ytitle("Parent religiosity effect on relationship") plotregion(color(white)) graphregion(color(white))  ///
saving(relate5, replace) ///
note(r =`corr') title("A. Religiosity effect on parenting and secular culture") legend(pos(6)order(1 "Linear Fit" 2 "Observed" )  region(style(none))  rows(1) span) 

corr B2  lngdp 
local corr : di %05.2g r(rho) 
twoway lfit (B2  lngdp) || scatter (B2  lngdp), mlabel(iso) ///
ytitle("Parent religiosity effect on relationship") plotregion(color(white)) graphregion(color(white))  ///
saving(relate6, replace) ///
note(r =`corr') title("B. Religiosity effect on parenting and GDP per capita") legend(pos(6)order(1 "Linear Fit" 2 "Observed" )  region(style(none))  rows(1) span) 

**FIGURE 3
grc1leg relate1.gph relate2.gph relate3.gph relate4.gph   ///
, imargin(2 2 2 2) row(3) col(2) legendfrom(relate1.gph) iscale(*.75) plotregion(color(white)) graphregion(color(white))  
graph export "Fig3_Predict_Relationship_Plot.png", width(1800) height(1200)  replace

**FIGURE 4
grc1leg relate5.gph relate6.gph   ///
, imargin(2 2 2 2) row(3) col(2) legendfrom(relate5.gph) iscale(*.75) plotregion(color(white)) graphregion(color(white))  
graph export "Fig4_Predict_Religion_Effect_Plot.png", replace width(1800) height(1200) 

*******
*TABLE--**Correlation STATS***
********
*install ci2 package if not installed
use "${clean_data}\gfs_cleaned_merged_database.dta", clear
merge 1:1 country_name using "${output}\religious_diversity_by_country.dta"

cor B1 B2 PCRQ wps_education wps_employment wps_fin_inclusion wps_cell_phone wps_parliament wps_no_legal_discrimination wps_justice wps_mortality_ratio wps_son_bias wps_partner_violence wps_safety wps_terrorism wps_conflict
cor B1 PCRQ  v2x_polyarchy zwvs_traditional zwps_index lngdp zmort_rate thriving HOPE HAPPY ACCEPTED_BY_GOD FINANCIAL parent_relig_index self_relig_index parents_married parent_died  flove_no_dad mlove_no_mom child_poverty was_abused
cor B2 PCRQ  v2x_polyarchy zwvs_traditional zwps_index lngdp zmort_rate thriving HOPE HAPPY ACCEPTED_BY_GOD FINANCIAL parent_relig_index self_relig_index parents_married parent_died  flove_no_dad mlove_no_mom child_poverty was_abused
cor DIV B1 B2 PCRQ self_relig_index zwvs_traditional

foreach x in thriving HOPE HAPPY v2x_polyarchy zwvs_traditional zwps_index lngdp ///
parent_relig_index self_relig_index B1 B2 T1 T2 {
ci2 PCRQ `x', corr
gen R1`x'=r(rho)
gen LB1`x'=r(lb)
gen UB1`x'=r(ub)
gen N1`x'=r(N)
}

foreach x in  PCRQ thriving HOPE HAPPY v2x_polyarchy zwvs_traditional zwps_index  ///
parent_relig_index self_relig_index B1 B2 T1 T2  {
ci2 lngdp `x', corr
gen R2`x'=r(rho)
gen LB2`x'=r(lb)
gen UB2`x'=r(ub)
gen N2`x'=r(N)
}

foreach x in PCRQ thriving  HAPPY v2x_polyarchy zwvs_traditional zwps_index lngdp ///
parent_relig_index self_relig_index B1 B2 T1 T2  {
ci2 HOPE `x', corr
gen R3`x'=r(rho)
gen LB3`x'=r(lb)
gen UB3`x'=r(ub)
gen N3`x'=r(N)
}

foreach x in PCRQ thriving HOPE HAPPY v2x_polyarchy zwvs_traditional zwps_index lngdp  ///
parent_relig_index self_relig_index  {
ci2 B1 `x', corr
gen R4`x'=r(rho)
gen LB4`x'=r(lb)
gen UB4`x'=r(ub)
gen N4`x'=r(N)
}

foreach x in PCRQ thriving HOPE HAPPY v2x_polyarchy zwvs_traditional zwps_index  lngdp ///
parent_relig_index self_relig_index  {
ci2 B2 `x', corr
gen R5`x'=r(rho)
gen LB5`x'=r(lb)
gen UB5`x'=r(ub)
gen N5`x'=r(N)
}

keep in 1
keep R1* R2* R3* R4* R5* LB* UB* 
gen n=1
reshape long R1 R2 R3 R4 R5 LB1 LB2 LB3 LB4 LB5 UB1 UB2 UB3 UB4 UB5, i(n) j(var) string
drop n
ren R1 COR_PCRQ
ren R2 COR_GDP
ren R3 COR_Flourishing
ren R4 COR_B1
ren R5 COR_B2

drop if var=="T1" 
drop if var=="T2"

export delimited "Country_Level_Correlations_CIs.csv", replace

